#if defined(OCL_KHR_FP64)
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#define MDS_F(X) X
#else
#define MDS_F(X) X##f
#endif


//==============================================================================
//
// Fortran anint in OpenCL...
//
//==============================================================================
Real anint( Real u )
{
	return floor( u + MDS_F(0.5) );
}

//==============================================================================
//
// PBC in OpenCL: -L/2 <= x <= L/2
//
//==============================================================================
Real pbc1( Real x, Real L )
{
	return x - L * anint( x/L );
}

//==============================================================================
//
// Vectorial PBC over xyz
//
//==============================================================================
Real4 pbc_xyz( const Real4 r, const Real4 box )
{
	return (Real4)( pbc1(r.x,box.x), pbc1(r.y,box.y), pbc1(r.z,box.z),  0 );
}



//==============================================================================
//
// Velocity Verlet data mapping
//
//==============================================================================
#define dt      (data.s0)
#define half_dt (data.s1)

//==============================================================================
//
// Velocity Verlet first part: must left all the .w fields untouched !
// 
//==============================================================================
__kernel void predict(__global   Real4 *position, 
					  __global   Real4 *velocity,
					  __global   Real4 *acceleration,
					  const      Real2  data,
					  const      Real4  box)
{
	const size_t    i   = get_global_id(0);
	__global Real4 *pos = position+i;
	__global Real4 *vel = velocity+i;
	Real4           r   = (Real4)((*pos).xyz,0);
	Real4           v   = (Real4)((*vel).xyz,0);
	Real4           a   = (Real4)(acceleration[i].xyz,0); //TODO: necessary ?
	
	v += half_dt * a;
	r += dt      * v;
	r = pbc_xyz(r,box);
	(*pos).xyz = r.xyz;
	(*vel).xyz = v.xyz;
	
}

//==============================================================================
//
// Velocity Verlet second part: must left all the .w fields untouched !
//
//==============================================================================
__kernel void correct(__global   Real4 *position, 
					  __global   Real4 *velocity,
					  __global   Real4 *acceleration,
					  const      Real2  data,
					  const      Real4  box)
{
	const size_t i      = get_global_id(0);
	__global Real4 *vel = velocity+i;
	// Real4           r   = (Real4)(position[i].xyz,0);
	Real4           v   = (Real4)((*vel).xyz,0);
	Real4           a   = (Real4)(acceleration[i].xyz,0); //TODO: necessary ?
	
	v += half_dt * a;
	(*vel).xyz = v.xyz;
}

//==============================================================================
//
// Lennard-Jones data mapping over a Real8
//
//==============================================================================
#define epsilon   (Vargs.s0)
#define sigma     (Vargs.s1)
#define sigma_sq  (Vargs.s2)
#define epsilon4  (Vargs.s3)
#define epsilon24 (Vargs.s4)

//==============================================================================
//
// Full Lennard Jones Dynamics
//
//==============================================================================
__kernel void dynamics(__global const  Real4 *position,
					   __global const  Real4 *velocity,
					   __global        Real4 *acceleration,
					   const           Real8  Vargs,
					   const           Real4  box,
					   __global        Real2 *output,
					   __local         Real4 *pcache,
					   __local         Real2 *vwcache
					   )
{
	const size_t gi      = get_global_id(0);   //-- global work_item index
	const size_t li      = get_local_id(0);    //-- local  work_item index
	const size_t gn      = get_global_size(0); //-- total size
	const size_t ln      = get_local_size(0);  //-- group size
	const size_t nblocks = gn / ln;            //-- number of groups
	
	const Real4  p   = position[gi];
	const Real4  r_i = (Real4)(p.xyz,0);
	Real4        f_i = (Real4)(0,0,0,0);
	Real         W   = 0; //! partial virial
	Real         V   = 0; //! partial potential
	
	//--------------------------------------------------------------------------
	//
	// Loop over all the work groups
	//
	//--------------------------------------------------------------------------
	for( size_t iblock=0; iblock < nblocks; ++iblock )
	{
		//----------------------------------------------------------------------
		//
		// cache #ln particle positions from the work group #iblock
		//
		//----------------------------------------------------------------------
		size_t start  = iblock * ln;
		pcache[li]    = (Real4)(position[start + li].xyz,0);
		barrier(CLK_LOCAL_MEM_FENCE);
		
		//----------------------------------------------------------------------
		//
		// loop over all particles in cache
		//
		//----------------------------------------------------------------------
		for( size_t j=0; j < ln; ++j )
		{
			const size_t  gj = start++; //! global index of particle j
			if( gj != gi )
			{
				Real4 r_ij = pcache[j] - r_i;
				Real4 f_ij;
				Real  V_ij;
				
				r_ij = pbc_xyz( r_ij, box );
				
				
				/** TODO: insert potential/foce code **/
				const Real  r2   = dot(r_ij,r_ij);
				const Real  q2   = sigma_sq/r2;
				const Real  q6   = q2 * q2 * q2;
				const Real  q12  = q6*q6;
				f_ij = ( epsilon24*( q12+q12 - q6)/r2)  * r_ij;
				//V_ij = epsilon4 * ( q12 - q6 );
				
				
				/** update force */
				f_i -= f_ij;
				//V   += V_ij;
				
				/** update virial **/
				//W += dot( r_ij, f_ij );
			}
			
		}
	}
	
	//--------------------------------------------------------------------------
	//
	// compute acceleration
	//
	//--------------------------------------------------------------------------
	f_i *= p.w;
	acceleration[gi].xyz = f_i.xyz;
	return;
	
	//--------------------------------------------------------------------------
	//
	// cache (V,W) for this work_item
	//
	//--------------------------------------------------------------------------
	vwcache[li] = (Real2)(V,W);
	barrier(CLK_LOCAL_MEM_FENCE);
	
	//--------------------------------------------------------------------------
	//
	// compute ONCE the partial (V,W) for the work_group
	//
	//--------------------------------------------------------------------------
	if( 0 == li )
	{
		//----------------------------------------------------------------------
		// local partial sum
		//----------------------------------------------------------------------
		Real2 vw = vwcache[0];
		for( size_t j=1; j < ln; ++j )
		vw += vwcache[j];
		
		//----------------------------------------------------------------------
		// global store
		//  gi/ln is the same for any work item in the work_group !
		//----------------------------------------------------------------------
		output[ gi/ln ] = vw;
	}
	
}

//==============================================================================
//
// a simple task to finish the (V,W) sum
//
//==============================================================================
__kernel void sum_VW( __global const Real2 *output, __global Real2 *VW, int n )
{
	Real2 vw = (Real2)(0,0);
	while( n-- > 0 )
	{
		vw += output[0];
		++output;
	}
	vw   *= MDS_F(0.5);   // since the terms are computed twice
	VW[0] = vw;           // done
}


